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ABSTRACT 

The Kelvin-Helmholtz instability is well known to be capable of converting well- 
ordered flows into more disordered, even turbulent, flows. As such it could represent 
a path by which the energy in, for example, bowshocks from stellar jets could be con- 
verted into turbulent energy thereby driving molecular cloud turbulence. We present 
the results of a suite of fully multifluid magnetohydrodynamic simulations of this in- 
stability using the HYDRA code. We investigate the behaviour of the instability in a 
Hall dominated and an ambipolar diffusion dominated plasma as might be expected 
in certain regions of accretion disks and molecular clouds respectively. 

We find that, while the linear growth rates of the instability are unaffected by 
multifluid effects, the non-linear behaviour is remarkably different with ambipolar 
diffusion removing large quantities of magnetic energy while the Hall effect, if strong 
enough, introduces a dynamo effect which leads to continuing strong growth of the 
magnetic field well into the non-linear regime and a lack of true saturation of the 
instability. 

Key words: mhd - instabilities - ISMxlouds - ISM:kinematics and dynamics 



1 INTRODUCTION 



The Kelvin-Helmholtz (KH) instability is an important in- 
• stability in almost any system involving fluids: it can occur 
I anywhere that has a velocity shear. In astrophysical plas- 
. mas the KH instability can provide the means of producing 
■ turbulence in a medium or the mixing of material between 

two boundary layers. 

The KH instability has been studi ed in a variety of 
' astrophysical systems, from solar winds ( Amerstorfer et al.l 



l2007 f: 'Bettarini et al."2006'; 'Hascgawa ct aLl bOO^ I and pul- 



sar w inds (Bucciantini & Del Zanna 2006) to thermal flares 
I Venter fc MeintiesI I200Q V Due to its ability to drive mix- 
ing and turbulence, the KH instability has been con- 
sidered relevant in protoplan etary disks (jJohansen et al.l 
120061 : iGo mez & Ostriker 2003), accretion disks and magne- 
tospheres (Li fc Naraya n 200^1), and other jets and outflows 
(jBaty fc Keppend 120061 ). Generally speaking, the assump- 
tions of ideal magnetohydrodynamics (MHD) have been 
used in order to simplify the system of equations to be 
solved. These assumptions are, however, not always valid. 
Weakly ionised plasmas, for example, contain a large frac- 
tion of neutral particles as well as a number of charged parti- 
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cle fluids with differing physical characteristics. Interactions 
between the various species can introduce non-ideal effects. 
Ambipolar dissipation and the Hall effect are two non-ideal 
effects that can greatly influence the development of the 
KH instability in a system by altering the dynamics of the 
plasma and the evolution of the magnetic fleld. Astrophys- 
ical examples of such weakly ionised systems in clude dense 
molecular clouds (e.g. ICiolek fc Robergel |2002|) and accre - 
tion disks around young stellar objects fe.g. IWardlelll999l '). 
In these systems, the relevant length scales are such that 
non-ideal effects can play an imp ortant role (|Wardlell2004al : 
iDownes fc O'Sullivan|[2009l . l201ll ) . 

Many authors have investigated the role of the KH 
instability in both magnetised and unmagnetised astro- 
ph ysical flow s (e.g. 'Frank et aLlll996l : iMala goU ct al lll996l: 
jHa rdee et al.l [l997; Downo s fc Ravi 1 19981 : Kcppcn s et al.l 
[l999). Most of these studies have investigated the KH in- 
stability in the context of either hydrodynamics or ideal 
MHD. We know that non-ideal effects are important in 
mol ecular clouds at len g th scales below about 0.2 pc (e.g. 
Ois hi fc Mac Low! l2006l : iDownes fc O'SuUivanl .2009) and 
hence it is of interest to explore the KH instability in the 
context of either non-ideal MHD or, preferably, fully multi- 
fluid MHD. In more recent years the empha sis of KH stud- 
ies has been on including non-ideal effects. iKeppens et al.l 
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l|l999l ) studied both the hnear growth and subsequent non- 
hnear saturation of the KH instability using resistive MHD 
numerical simulations. The inclusion of diffusion allowed for 
magnetic reconnection and non-ideal effects were observed 
through tearing instabilities and the formation of magnetic 
islands. 

IPalotti et al.l (|2008l ) also carried out a series of simu- 
lations using resistive MHD. They found that, following its 
initial growth, the KH instability decays at a rate that de- 
creases with decreasing plasma resistivity, at least within the 
range of of resistivities accessible to their simulations. They 
also found that magnetisation increases the efficiency of mo- 
mentum transport, a nd that the transport increases with 
decreasing resistivity. iBirk fc WiechenI (|2002l ) examined the 
case of a partially ionised dusty plasma, using a multifluid 
approach in which collisions could be included or ignored. 
They found that collisions between the neutral fluid and 
dust particles could lead to the stabilisation of KH modes 
of particular wavelengths. The unstable modes led to a sig- 
nificant local amplification of the magnetic field strength 
through the formation of vortices and current sheets. In the 
nonlinear regime they observed the magnetic flux being re- 
distributed by magnetic reconnection. It was suggested that 
this could be applicable to dense molecular clouds and have 
important implications for th e magnetic flux loss problem 
l|Umebavashi fc Na kano"l990Y 

A comprehensive study was carried out by IWiechenI 

1(2003) which demonstrated the effect of dealing with the 
plasma using a multifluid scheme. This study focused on 
the effect of varying the properties of the dust grains. The 
results of the simulations led to the conclusions that more 
massive dust grains have a stabilising effect on the system 
while higher charged numbers have a destabilizing effect. It 
was found that there is no significant dependence on the 
charge polarity of the dust. 

In a linear study of stellar outflows I Wat son et al.l (|2004l ) 
described how the charged and neutral fluids are affected 
differently by the presence of a magnetic fleld. This study 
is carried out using parameters chosen to reflect those of 
molecular clouds, and so is particularly relevant to our own 
study. The principal result of this paper is that for much 
of the relevant parameter space, neutrals and ions are suf- 
ficiently decoupled that the neutrals are unstable while the 
ions are held in place by the magnetic field. Since the mag- 
netic field is frozen to the ionised plasma, it is not tangled 
by the turbulence in the boundary layer. The authors pre- 
dict that with well-resolved observations, there should be a 
detectably narrower line profile in ionised species tracing the 
stellar outflow compared with neutral species, since ionised 
species are not participating in the turbulent interface with 
the ambient interstellar medium. The paper also includes a 
study of the growth rate of the instability. It is found that 
at short length scales, the growth rate is well approximated 
by the growth rate of the hydrodynamic system. At larger 
scales and for super- Alfvenic flows, the fastest growing mode 
is e qual to that of the idea l MHD case. 

IShadmehri fc Downea (|2008l ) carried out an analytical 
study of the Kelvin-Helmholtz instability in dusty and par- 
tially ionised outflows. They investigated primarily the ef- 
fect of the presence of dust particles by varying their mass, 
charge and charge polarity. It was found that as the charge of 
the grain increased, the growth timescales also increased, im- 



plying a stabilising effect on the system. The stability of the 
system was also examined for dependence on the mass of the 
dust particles. It was found that for stronger magnetic fields, 
this did not affect the stability of the system. However, for 
weaker magnetic fields, the larger dust particles had a sta- 
bilising effect on the growing modes. T his was in agree ment 
with previous laborato ry experiments (jLuo et al.ll200ll ) and 
numerical simulations ( Wiec henll200q ) . Finally, as the mag- 
netic field strength increased, the growth timescale of the 
unstable modes at a particular perturbation wavelength de- 
creased. By examining the combinations of the wavelength of 
the perturbation u sed, and the resultant growth timescales 
of the instability, IShadmehri fc Downed 1 20081 ') concluded 
that the Kelvin-Helmholtz instability is a possible candidate 
for causing the formation of some of the physical structures 
observed in molecular outfiows from young stars. 

In this paper we perform numerical simulations of the 
complete evolution of the Kelvin-Helmholtz instability in a 
weakly ionised, multifluid plasma including both its linear 
development, saturation and its subsequent behaviour. We 
include the physics of the Hall effect, ambipolar diffusion 
and parallel resistivity in these simulations and we analyse 
the role each of the former two effects on the development 
of the instability. 

The aim of this work is to investigate the influence of 
multifluid effects on the growth and saturation of the KH 
instability. In order to develop a full understanding of the 
roles of the various non-ideal effects, in particular the Hall 
effect and ambipolar diffusion, we run simulations with pa- 
rameters chosen to simulate very high, medium and very 
low magnetic Reynolds number systems and with parame- 
ters chosen to ensure ambipolar-dominated fiows and Hall- 
dominated flows. We focus on gaining an understanding of 
the general characteristics of the KH instability in weakly 
ionised flows. A detailed study of this instability in the spe- 
cific context of molecular clouds, and which is of interest 
from the point of view of turbulence generation by stellar 
outfiows, is the subject of a future work. 

In section[2]we outline the numerical and physical model 
employed, in sectionOwe discuss how we analyse the results 
of our simulations while in section [4] we detail the results in 
both the linear and non-linear regimes, separating out the 
effects of ambipolar diffusion and the Hall effect in order to 
more fully understand the infiuence of each. 



2 NUMERICAL APPROACH 

The simulations describ ed in this work are performed us - 
ing the HYDRA code (|0'Sullivan fc Downed I2OO6I . |2007| ) 
for nmltifiuid magnetohydrodynamics in the weakly ionised 
regime. We further assume the flow is isothermal. The as- 
sumption of weak ionisation allows us to ignore the inertia 
of the charged species and allows us to derive a (relatively) 
straightforward generalised Ohm's law. The resulting sys- 
tem of equations, given below, incorporates finite parallel, 
Hall and Pederson conductivity and is valid in, for example, 
molecular clouds. In such regions the viscous lengthscales 
are much smaller than those over which nonideal effects are 
important. This leads to high Prandtl numbers and plasma 
flows in these regions can be considered to be effectively 
inviscid. In this work we examine low Mach number flows 
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which, taken in concert with the isothermal assumption, 
means that features in the flow such as shocks are unlikely 
to create regions of high ionisation. 

2.1 Multifiuid equations 

The HYDRA code solves the following equations for a sys- 
tem of A^ fluids. The simulations described in this paper 
consist of three fluids, indexed by i = for the neutral fluid 
and 1 = 1 and i = 2 for the electron and ion fluids respec- 
tively. The equations to be solved are 



^+V-(p,qO=0 (0<J^iV-l), 



^po^o 



dt 

dB 

dt 



+ V ■ (poqoqo + a^pol) = J X B, 

+ V ■ (qoB-Bqo) = -V x E', 



(1) 
(2) 
(3) 



aiPi{E + ciiXB)+pipoKio{(io-cii) = (1 sj i sj iV-l),(4) 
V • B = 0, (5) 

VxB = J, (6) 



ciipi = 0, 



(7) 



N-1 

E 

i=l 

where pi, q^, B, and J are the mass densities, velocities, 
magnetic field and current density, respectively, a denotes 
the sound speed, and ai and Kio are the charge-to-mass ra- 
tios and the collision coefficients between the charged species 
and the neutral fluid, respectively. 

These equations lead to an expression for the electric 
field in the frame of the fluid, E', given by the generalised 
Ohm's Law 

E' = Eo+Eh + Ea, (8) 

where the components of the field are given by 

Eo = (J • ao)ao, (9) 

Eh = J X an, (10) 

Ea = -(J X an) X an, (11) 

using the definitions ao = /oB, an = /hB, aA = /aB, 
where /o = y^/B, fu = th/B and /a = y^/B. The 
resistivities given here are the Ohmic, Hall and ambipolar 
resistivities, respectively, and are defined by 



ro 



o-Q 



o-H 





ai + al' 










OK 










2 1 2 ' 

0-^ + crj^ 




wher< 


e the conductivities 


are 


given 


by 


o-o = 


^ N-1 
■ - ^ Q,p,ft, 
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0"H = 


N-1 
1 \ ^ OLipi 

B 2^ l + Pf' 

i—l 









(12) 
(13) 
(14) 

(15) 
(16) 



JV-l 

1 v^ OiipiPi 

i—l 



(17) 



where the Hall parameter fii for a charged species is given 
by 

a^B 



ft = 



(18) 



Kiopo 

To solve these equations numerically we use three dif- 
ferent operators: 

(i) solve equations (1), (2), (3), including the restriction 
of equation (5) and for i = 0, using a standard second or- 
der, finite volume shock-capturing scheme. Note that for this 
operator the resistivity terms in equation (3) are not incor- 
porated. Equation (5) is incorporated using the method of 
Dedner |E)edner et al.ll2002l ). 

(ii) Incorporate the resistive effects in equation (3) us- 
ing super-time-stepping to accelerate the ambipolar diffu- 
sion term and the Hall Diffusion Scheme to deal with the 
Hall term. 

(iii) Solve equations (4) for the charged species velocities 
and use these to update equation (1) (with i = 1, . . . ,N — 1) 

These operators are applied using Strang operator splitting 
in order to maintain the second order accuracy of the overal l 
schem e. We refer the reader to lO'Sullivan fc Downed (|2006l . 
[20071) for a more detailed description. 

2.2 Initial conditions 

The simulations are carried out on a 2.5 D slab grid in the 
xy-plane. The grid consists of 6400 x 200 x 1 cells, in the x, 
y, and z directions respectively. This resolution was chosen 
on the basis that it reproduces t he initi al linear growth of 
the ideal MHD system in Kcppcn s et al.l (11999). Resolution 
studies were performed to confirm the resolution as being 
appropriate (see Sect. I4.2~l1 and l4.3.1|) . 

The initial set-up used was that of two plasmas fiowing 
anti-parallel side-by-side on a grid of size x — [0, 32L] and 
y — [0, L]. The plasma velocities are given by -|-^ and — ^ 
in the y-direction, with a tangential shear layer of width 2a 
at the interface at x — 16L. This velocity profile is described 
by 



vo ^ 



— tanh I — 



16L\ , 

— y- 

a I 



(19) 



The width of the shear layer is chosen to be j- = 0.05, or 
approximately 20 grid zones. The magnetic field is initially 
set to be uniform and aligned with the plasma flow. 

The initial background for all three fluids in the system 
is now an exact equilibrium. The initial neutral velocity fleld, 
Vb is then augmented with a perturbation given by 



&V:c — SVq sin{—kyy) exp 



(x- 16Lf 



(20) 



where dVo is set to 10 * Vq. The wavelength of the per- 
turbation is set equal to the characteristic length scale. 



A 



L, so that a single wavelength flts exactly into 



the computational domain. This maximises the possibility 
of resolving structures t hat are small relat ive to the initial 
perturbed wavelength, (jFrank et al.l 119961 ). The perturba- 
tion attenuation scale is chosen so that it is larger than the 
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Fluid 


Density 


Kio 


Oi 


1 

2 


2.84 X 10"^^ 
2 X 10"'^ 


2 X lO" 
1.42 X 10^" 


-1 X lO^'^ 
1.42 X 10" 



Table 1. The density, collision coefficient (Kio) and charge-to- 
mass ratios (a^ ) for each of the charged fluids in simulation fuU- 
low-hr (see table [2]|. These parameters are modified to vary the 
resistivities as necessary for the other simulations in this work. 
See text. 



shear layer, but small enough so that the instability can be 
assu med to interact onl y minimally with the a;-boundaries 
2008|'). and is set using 



(see Palotti et al.l 



iKeppens et al.ll 19991 : iPalotti et al. 2008). 



f = 0.2 



(see 



The physical parameters are then chosen using nor- 
malised, dimensionless quantities. The wavenumber ky is 
chosen to b e 27r in order to max imise the growth rate of the 
instability (JKeppens et al.ll 19991 ). This normalises the length 
scale of the simulation so that L = 1. The timescale is then 
normalised by the sound speed, so that Cb = -^ = 1. The 
mass scale is chosen such that the initial mass density is set 
to unity, po = 1- In the isothermal case the adiabatic index 
7 = 1, and this gives us a sound speed Cs = ^JTSK — \^ 
and the initial pressure is therefore also equal to unity, 
po = 1. A transonic flow is chosen with sonic Mach number 
Ms = ^ = 1, so that the plasma has velocity ±^ ~ ±0.5. 
In order for the KH insta bility to grow, it is re quired that 
the flow be super- Alfvenic (|Chandrasekhaiill96l| ). and so for 
this study, an Alfven Mach number M\ = — 2- = 10 is cho- 

sen. This sets the Alfven velocity, va = \ —^ ~ 0.1, and 
the initial magnetic field strength to Bo — 0.1. 

We use periodic boundary conditions at the high and 
low y boundaries. Since we wish to study not only the initial 
growth phase of the instability, but also its subsequent non- 
linear behaviour we must ensure that waves interacting with 
the high and low x boundaries do not reflect back into the 
domain to influence the dynamics. Several test simulations 
for various parameters have shown that a large width of 32 
is necessary to ensure this. We use gradient zero boundary 
conditions at the high and low x and z boundaries. 

Finally we must choose the parameters describing the 
properties of our charged fluids in our multifluid s ystem. 
Our b asic parameter set is contained in table [T] See Ijoned 
l|2011f) for more details. 

Table [2] contains the nomenclature we will use for the 
rest of this paper when referring to the simulations. Each 
simulation is denoted by xxxx-yyyy-zz where the first set 
of characters denote the dominant resistivity (Hall or am- 
bipolar or "full" if both resistivities have the same mag- 
nitudes), the second set denote the level of the resistivity 
(low, medium or high) and the final two characters denote 
the resolution (low, medium or high resolution denoted by 
Ir, mr and hr respectively). The high resolution simulations 
(6400 X 200) took of the order of 1500 - 3000 core hours on 
a quad-core Xeon E5430 based system. 

Note that simulation fuU-low-hr is effectively an ideal 
MHD simulation and produced results virtually identical 
to mhd-zero-hr which is a true MHD simulation (see Sect. 
I4.1|) and which was used for comparison with previous liter- 
ature. 



Simulation 


Resolution 




ru 






VA 




mhd-zero-hr 


400 X 200 
















hd-zero-hr 


400 X 200 
















ambi-high-lr 


1600 X 50 


3.52 


X 10" 


-6 


3.5 


X 10" 


-2 


ambi-high-mr 


3200 X 100 


3.52 


X 10" 


-6 


3.5 


X 10" 


-2 


ambi-high-hr 


6400 X 200 


3.52 


X 10" 


-6 


3.5 


X 10" 


-2 


fuU-low-hr 


6400 X 200 


3.52 


X 10" 


-B 


3.5 


X 10" 


-B 


ambi-med-hr 


6400 X 200 


3.52 


X 10" 


-6 


3.5 


X 10" 


-3 


hall-high-lr 


1600 X 50 


3.52 


X 10" 


-2 


3.5 


X 10" 


-B 


hall-high-mr 


3200 X 100 


3.52 


X 10" 


-a 


3.5 


X 10" 


-B 


hall-high-hr 


6400 X 200 


3.52 


X 10" 


-2 


3.5 


X 10" 


-B 


hall-med-hr 


6400 X 200 


3.52 


X 10" 


-3 


3.5 


X 10" 


-B 



Table 2. The resolutions and (initial) resistivities of each of the 
simulations presented in this work. The Hall resistivity is denoted 
by rjj while the ambipolar resistivity is denoted by r a . Note that 
the Pederson conductivity is always very large and hence its as- 
sociated resistivity is always several orders of magnitude smaller 
than both th and ta. 



3 ANALYSIS 

In order to study the growth of the instability, the evolution 
of a number of parameters can be measured with time. In 
particular, we measure the transverse kinetic energy 



Ev 



and the magnetic energy 



— pv^ dx dy 



Eb = l f f {[B'.+B^y+ B,'] - Bo'} dxdy 

in the system where Bo is the magnitude of the magnetic 
field at t = 0. As the entire plasma flow is initially in the 
y-direction, with only a very small perturbation in the x- 
direction, any growth of Bk.i: is due to the growth of the 
instability. It is possible to determine the growth rate of the 
instability dir ectly from the growth of the transverse kinetic 
energy (Kep pens et al.lll99& l: the transverse kinetic can be 
expressed as 



Bk,i: = {po + 5p)(\vx\o + Sv^ 

= {po + Sp){5vx)' 

~ PoSv^ 

oc exp[2i(k • r — tot)] 



(21) 
(22) 
(23) 
(24) 



presuming the initial perturbation is proportional to exp[i(k- 
r—ujt)] . Hence the kinetic energy grows at a rate of 2ijj, where 
uj is the growth rate of the instability. 



4 RESULTS 

We start by describing the validation of the set-up used by 
comparing an ideal MHD simulation run using HYDRA with 
previously published literature. We then go on to discuss the 
behaviour of the KH instability in ambipolar-dominated and 
Hall-dominated flows respectively. 
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Figure 1. Plot of the log of the transverse kinetic energy with 
time for mhd-zero-hr. A linear function (dashed line) with slope 
2.63 has been over-plotted. 



Figure 3. Plot of the evolution of the transverse kinetic energy 
against time in the case of high ambipolar resistivity for ambi- 
high-lr, ambi-high-mr and ambi-high-hr. 




Figure 2. Plot of the evolution of magnetic energy with time for 
simulation mhd-zero-hr. 



4.1 Validation 

In order first to validate our set-up we examine our mhd- 
zero-hr simulation with HYDRA and determine the growth 
rate using the kinetic energy of motions in the x direction 
as described in Sect. (3] Figure [T] contains a plot of the log of 
the transverse kinetic energy as a function of time. At early 
times this growth is clearly exponential and can be fitted 
with a line of slope 2.63 implying a growth rate, normalised 
by the width of the shear layer and the initial relative veloc- 
ity, for the dominant mode of the KH instability of 0.1315. 
We can compare this with the value of t he growth rate calcu- 
lated analytically bv lMiura fc PritchettI ([1982) (their Figure 
4) at this wavenumber of 0.13. While comparisons between 
linear studies of incompressible flows, and numerical stud- 
ies of compressible flows are bound to differ to some extent, 
these results are seen to agree exceptionally well. 

We wish to examine not only the linear regime but also 
the non-linear regime. We compa re our results for the gr owth 
of magnetic energy with those of lMalagoli et al.1 (| 19961 ) (the 
upper panel of their Figure 5). Figure[2]contains a plot of the 
magnetic energy, calculated as J J i(_B^-f B^-f _B^j dxdy, as 
a function of time. The maximum magnetic energy reac hed 
in our simulations matches that of lMalagoli et al.l ()l996l ) to 
within 10%. Our simulation reaches saturation at a later 
time but the exact time of saturation depends on the initial 
amplitude of the perturbation and so this is not a concern. 

We are therefore confident of the behaviour of HYDRA 
in simulating the KH instability. We now move on to investi- 
gating the influence of multifluid MHD effects on the growth, 
saturation and non-linear behaviour of this instability. 



4.2 Ambipolar dominated flows 

We begin our study of the nmltifluid KH instability by 
choosing our fluid parameters to ensure our ambipolar resis- 
tivity is dynamically significant while minimising the Hall 
resistivity. This allows us to isolate the infiuence of the am- 
bipolar resistivity on the instability. In order to increase the 
ambipolar resistivity we change the value of the collision 
coefficient for species 2 so that K2,o is decreased by 3 or- 
ders of magnitude from that given in table [T] for simulation 
ambi-med-hr and by 4 orders of magnitude for ambi-high-hr. 
These alterations of ii'2,0 give values of ta of 3.5 x 10""^ and 
3.5 X 10~^ respectively, and (ambipolar) magnetic Reynold's 
numbers. Rem, of 2.84 x 10^ and 2.84 x 10^ respectively. 
These simulations are examined in comparison to the fuU- 
low-hr simulation. With a formal magnetic Reynolds num- 
ber 2.84 X 10^, the diffusion in this set-up is predominantly 
numerical, and as such, it is effectively an ideal MHD simu- 
lation. 



4-2.1 Resolution study 

In non-ideal MHD we must ensure that the length scales 
over which the diffusion of the magnetic field (or the whistler 
waves in the case of Hall dominated flows) must be resolved 
in order to properly track the dynamics of the system. To 
this end we perform a resolution study using simulations 
ambi-high-lr, ambi-high-mr and ambi-high-hr (see Table [2]). 

Figures |3] and |4] contain plots of the evolution of -Ek,i^ 
and _Eb for each of the simulations in our resolution study. It 
can be seen that the linear growth in ambi-high-lr is signif- 
icantly lower than the two other simulations. However, the 
linear behaviour is almost identical for ambi-high-mr and 
ambi-high-hr. The subsequent non-linear behaviour is simi- 
lar with only relatively small variations after i ~ 11. 

The results of this study indicate that a resolution of 
6400 X 200 is sufficient to capture the initial growth and 
saturation of the instability. Subsequently, the dynamics is 
captured at least qualitatively. 



4-2.2 The linear regime 

Generally speaking, the evolution of the KH instability in 
ideal MHD leads to a wind-up of both the plasma and the 
magnetic field at the interface between the two fluids, re- 
sulting in the "Kelvin's cat's eye" vortex. Multifluid effects 
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CD 



m 



0.0012 
0.0010 
0.0008 
0.0006 
0.0004 
0.0002 
0.0000 








6400x200x1. 
3200x100x1. 
1600x50x1 . 



20 



Figure 4. As figure [3l but for the perturbed magnetic energy. 





Figure 5. Plots of the magnitude and vector field of the magnetic 
field for fuU-low-hr, at time t = 8ts (upper panel) and for ambi- 
high-hr (lower panel). The magnetic field is wound up in a similar 
fashion to the velocity in the upper panel, but not the lower panel. 



alter the nature of this significantly. Figure [5] contains plots 
of the magnetic field a,t t = St^ which is the time at which 
the instability saturates for both fuU-low-hr (effectively ideal 
MHD) and ambi-high-hr. The difference in the morphology 
is clear. 

Given these striking differences at saturation it is in- 
teresting to investigate whether the linear growth rate is 
influenced by the addition of ambipolar diffusion. Figure [S] 
contains plots of the evolution of -Ek,i! with time for the 
various simulations. The linear growth rate remains almost 
unchanged with the addition of ambipolar diffusion. On 
the other hand, flgure [7] contains plots of Et as a func- 
tion of time. It is clear that the perturbed magnetic energy 
is strongly influenced, even well before saturation, by the 
presence of ambipolar diffusion. We will discuss this in more 
detail in section [4.2.31 




Figure 6. Plots of the log of the transverse kinetic energy 
for varying ambipolar resistivity. The thicker lines represent the 
models with higher ambipolar resistivity. The dashed line repre- 
sents hd-zero-hr. There is very little difference between the linear 
growths for the three magnetised cases. 
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Figure 7. The perturbed magnetic energy in the system is plot- 
ted against time in the three cases of varying ambipolar resistivity. 
The thicker lines represent the models with higher ambipolar re- 
sistivity. Higher ambipolar resistivity causes the magnetic field to 
experience less amplification, through diffusion and decoupling of 
the bulk fluid from the magnetic field. 



4-2.3 The nonlinear regime 

It is clear from figure [7] that the magnetic energy is strongly 
influenced by the presence of ambipolar diffusion. This is 
not too surprising as ambipolar diffusion, being a genuinely 
diffusive process (unlike the Hall effect) allows the magnetic 
field to diffuse relative to the bulk flow. 

As the collision rate between the ion and neutral fluids 
is decreased, the ion fluid decouples from the bulk fluid, and 
thus the magnetic fleld becomes decoupled from the bulk 
flow: the frozen-in approximation of ideal MHD is broken. 
As a result, the magnetic field is able to diffuse through the 
bulk fluid rather than being tied to it. In figure [71 diffusion 
can be identified as the cause of the decrease of the am- 
plification of the magnetic energy with time for increasing 
ambipolar resistivity. A cursory examination of the topology 
of the magnetic field for the low resistivity case (fuU-low-hr 
in figure[5} demonstrates clearly that there are regions in the 
flow which will be susceptible to diffusion: regions in which 
the magnetic field lines have been compressed and amplified. 

It is clear, again from figure [3 that there is a signifi- 
cant decrease in the growth of magnetic energy as a result 
of diffusion for even a moderate amount of ambipolar re- 
sistivity. The diffusion has an influence on the dynamics of 
the neutral fluid: when the magnetic energy has not been 



© 



RAS, MNRAS 000.[THT2l 



The Kelvin- Helmholtz instability in weakly ionised plasmas 7 





2.0-1 0-' 




1.5'10' 


> 
Q. 


1.0-10^ 


cy 


S.O'IO'" 









0.010 
0.008 
V 0.006 
E^ 0.004 
0.002 
0.000 



— 6400x200x1 
- 3200x100x1 
■-1600x50x1 








20 



Figure 8. Plots of the transverse kinetic energy of the ion fluid 
as a function of time for each of full-low-hr, ambi-med-hr, ainbi- 
high-hr and hd-zero-hr. The thicker hnes represent the models 
with higher ambipolar resistivity. The dashed line represents the 
hydrodynamic case. 



as strongly amplified, the field can no longer exert the same 
effect on the neutral and ion fluid. Examination of figure [6] 
reveals that the peak reached by E]^,^ increases with increas- 
ing ambipolar resistivity. The value reached tends toward 
the hydrodynamic limit for two reasons: increasing resistiv- 
ity implies increasing diffusion and hence a decreasing field 
strength, and it also implies less coupling between the mag- 
netic field and the neutral fiuid with increasing resistivity. 
As expected, the weaker magnetic field strength allows the 
vortex to become more rolled up (jFagancllo ct al. 2009). 

It can be shown that for the ambi-med-hr simulation, 
the increase in the wind-up of the bulk fiuid is due solely 
to the first source: the magnetic diffusion. The slight de- 
coupling of the ion and neutral fiuid is sufficiently high to 
allow magnetic diffusion, while still being sufficiently low to 
force the charged fiuids to behave in a manner similar to 
the bulk fiuid. This can be seen in a plot of the transverse 
kinetic energy of the ion fiuid, as shown in figure [8l With 
moderate amounts of ambipolar diffusion, the transverse ki- 
netic energy of the ion fluid (and electron fiuid) reaches a 
higher maximum, meaning that the fiuid experiences more 
wind-up. 

However, as the ambipolar resistivity is increased fur- 
ther, the ion fiuid becomes more decoupled from the neutral 
fiuid, and is instead infiuenced primarily by the dynamics of 
the magnetic field. This significant decoupling occurs only 
for magnetic Reynolds number lower than Rem ~ 100. As 
the magnetic field therefore undergoes less wind-up than be- 
fore, so does the ion (and electron) fiuid. This decrease in the 
winding-up of the charged fiuid velocity field is the cause of 
the decrease in the transverse kinetic energy of the ion fiuid, 
as is seen in ambi-high-hr in figure O 



4.3 Hall dominated flows 

We now turn our attention to the likely infiuence of the 
Hall effect on the KH instability in multifiuid MHD fiows. 
In order to attain the resistivities we want (see table ^ 
we reduce the charge-to-mass ratio of the ion fiuid, causing 
its Hall parameter to become smaller, causing the ion fiuid 
dynamics to more closely resemble the bulk (neutral) fiuid 
dynamics. The electrons are, however, still well-tied to the 
field lines and so a relative drift emerges between the ion 



Figure 9. Plot of the transverse kinetic energy against time for 
simulations hall-high-lr, hall-high-mr and hall-high-hr. 



and electron fiuids leading to a current perpendicular to the 
magnetic field and hence the Hall effect. 



4-3.1 Resolution study 

As in the ambipolar resistivity study, a resolution study is 
carried out to ensure that the small-scale non-ideal dynam- 
ics are captured. For this purpose, simulations are again run 
at three different resolutions (see Table [U . These simula- 
tions are run with the highest level of Hall resistivity to 
ensure that the smallest-scale dispersive effects are in place 
when examining whether they are sufficiently resolved. The 
inclusion of the Hall term in the dispersion relation in prin- 
ciple allows for the introduction of waves with a signal speed 
which tends towards infinity as their wavelength tends to- 
wards zero. While the Hall term is handled in the equa- 
tions by the H YDRA code using the explic i t Hal l Diffusion 
Scheme (HDS) (|0'Sullivan fc Downes|[200^ . l2007l ). the code 
naturally does not resolve these waves of vanishing wave- 
length. 

As has been seen, the introduction of non-ideal effects 
does not tend to greatly infiuence the linear growth rate 
of the instability. As a result, the growth rate does not pro- 
vide a good means of measuring convergence with increasing 
resolution. The nonlinear evolution of the transverse kinetic 
energy, _Bk,a; is, however, strongly infiuenced by the non-ideal 
effects. We do not examine the evolution of the perturbed 
magnetic energy as, in the case of high Hall resistivity, it no 
longer demonstrates a simple growth to an initial maximum. 
The evolution of the transverse kinetic energy for each sim- 
ulation is plotted in figure [21 It can clearly be seen that the 
simulations have started to converge at higher resolutions. 
While there is a notable difference between the two simula- 
tions of lower resolution, the gap closes significantly in the 
comparison between the two simulations of higher resolu- 
tion. In particular, the initial maxima of _Ek,a; are almost 
identical in simulations hall-high-mr and hall-high-hr. 

We are, therefore, confident that the dynamics in the 
hall-hr are well resolved and that our conclusions as to the 
physical processes occurring are well-founded. 



4-3.2 The linear regime 

Figures fTUl and [TT] contain plots of -Ek,!! and Et as a function 
of time for simulations hall-med-hr and full-low-hr. In the 
linear regime, the infiuence of the Hall effect on the growth 
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Figure 10. The log of the transverse kinetic energy in the sys- 
tem is plotted against time for the fuU-low-hr (thin line) and 
hall-med-hr (thick line). The dashed line represents hd-zero-hr. 
There is very little difference between the linear growths for the 
magnetised systems. 




Figure 11. The perturbed magnetic energy in the system for 
both full-low-hr (thin line) and hall-med-hr (thick line) are plot- 
ted against time. No significant difference can be seen between 
the maxima reached in each case. 



rate of the instability, and on the magnetic and kinetic en- 
ergy in the system at saturation can be seen to be negligible. 
Hence we can conclude that neither ambipolar diflfusion nor 
the Hall effect influence the energetics of the system in the 
linear regime. 

We expect the Hall effect to re-orient the magnetic field 
out of the xy -plane in which it resides at f = 0. Figure [T5] 
contains plots of Eb, the perturbed magnetic energy in the 
xy-p\a.ne and the perturbed magnetic energy in the z direc- 
tion as functions of time. It is clear that there is some growth 
of the magnetic field in the z direction. Interestingly, at satu- 
ration the magnetic energy in the a;j/-plane is noticeably less 
than El, and yet the overall value of Ei, is the same as that 
derived from the full-low-hr (i.e. quasi-ideal MHD) simula- 
tion. An interesting point to note abou t this is that, whereas 
it has been shown (jJones et al.lll997i ) that the strength of 
the field in the direction of the initial flow is what is impor- 
tant in determining the effect of the magnetic field on the 
KH instability, here we can see that the strength of the field 
perpendicular to this plane appears to play a role also. 



4-3.3 The nonlinear regime 

Following the initial linear growth of the instability, the sys- 
tem experiences a period of transferring energy back and 
forth between the magnetic field and velocity field, as in 



Figure 12. The total perturbed magnetic energy (solid line) in 
the system is plotted against time for hall-med-hr. Also plotted 
for this simulation are the evolutions of the perturbed magnetic 
energies in the xy-pla,ne only (dashed line) and in the ^-direction 
(dot-dashed line) . While the maximum reached by the total mag- 
netic energy is the same as in full-low-hr, the magnetic field has 
experienced some re-orientation into the ^-direction. 




Figure 13. Plot of the total perturbed magnetic energy (dot- 
dashed lines) and transverse kinetic energy (solid lines) with time 
for the full-low-hr (thin lines) and hall-med-hr (thick lines) sim- 
ulations. 



full-low-hr (see figure [TSj. Interestingly, the amplitude of 
the oscillations - i.e. the amount of energy being transferred 
between motion and the magnetic field - is larger in hall- 
med-hr than in full-low-hr. This may be due to re-orientation 
of the magnetic field out of the plane of the instability and 
hence a reduction in (numerical) reconnection. It is worth re- 
calling here that the Hall effect, although it appears similar 
to a diffusion term in the induction equation, is a dispersive 
effect which conserves magnetic energy. 

Perhaps the most important difference between full-low- 
hr and hall-med-hr is that while the kinetic energy in the 
y-direction gradually tends towards a constant value in the 
low resistivity case, it continues to steadily decay in hall- 
med-hr. The conclusion is that in full-low-hr, the instability 
has completed its growth and is returning to a quasi-steady 
state. In hall-med-hr however, the instability is continuing 
to consume the parallel kinetic energy available to it and 
the instability undergoes a further stage of development as 
a result of the inclusion of the Hall effect. 

Since the magnetic field gains a component in the z 
direction in hall-med-hr due to the Hall effect it is interest- 
ing to examine the behaviour of the kinetic energy in the 
z direction also. Figure [TJ] contains plots of the transverse 
kinetic energy and the kinetic energy in the z direction (i.e. 
_Ek,z = \pv1). Clearly the transverse kinetic energy grows 
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Figure 14. Plot of the transverse kinetic energy (solid line) and 
kinetic energy in the z-direction (dashed line) with time for hall- 
med-hr. It can be seen that at later times, the growth of energy 
in the 2-direction has become significant. 



rapidly during the linear development of the instability but 
i?k,z also grows and, eventually, even becomes larger than 
_Ek,i:. Figure [151 contains plots of the total energy (kinetic 
and magnetic) in fuU-low-hr, hall-med-hr and hall-high-hr. 
Somewhat surprisingly we find that the hall-med-hr simula- 
tion loses energy somewhat faster than fuU-low-hr. This is, 
on the face of it, a little puzzling since the Hall effect does 
not itself dissipate magnetic energy. If we consider the hall- 
high-hr simulation we find that the total energy is roughly 
constant - it behaves roughly the same as hd-zero-hr. In 
fact, the magnetic energy in the hall-high-hr simulation be- 
haves qualitatively differently to the hall-med-hr. Figure [161 
contains plots of the perturbed magnetic energy as a func- 
tion of time for hall-high-hr, and fuU-low-hr. It is clear that 
something dramatic is happening in the hall-high-hr simula- 
tion. To gain insight into this, in figure [TTl we plot the total 
perturbed magnetic energy and the magnetic energy in the 
a;2/-plane and in the z-direction for simulation hall-high-hr. 
The growth of the magnetic field in the z direction is signifi- 
cant, as we expect from a system with the Hall effect. In fact, 
the growth of the magnetic field in the sy-plane is faster in 
the hall-med-hr case. This isn't too surprising as, in order to 
increase the Hall effect, the coupling between the electrons 
and the neutrals is much weaker than that between the ions 
and neutrals. Hence, while the neutrals and ions generate 
the usual cat's eye vortex, the electrons (and the magnetic 
field) do not. Figure [191 demonstrates these morphological 
differences between the various fluids. 

Furthermore, in the case of fuU-low-hr, the magnetic 
field in the xy-]Aa,ne grows to such an extent that it op- 
poses further wind-up of the fluid in the vortex. In the hall- 
high-hr case this does not happen as magnetic energy is re- 
distributed to the 2 component of the field which does not 
oppose this wind-up. Hence the vortex is not destroyed in the 
non-linear regime in the hall-high-hr case. This is the main 
difference between the hall-med-hr and hall-high-hr simula- 
tions. In the non-linear regime, then, the ion fluid remains 
spinning in the KH vortex while the electron fluid remains 
tied to the magnetic field. This maintains a velocity differ- 
ence between the electrons and ions which, in turn, causes 
the magnetic energy to be further re-distributed to the z 
component. In this way the Hall effect, if strong enough, in- 
troduces strong dynamo action into the KH instability. This 
dynamo behaviour, which is not possible in a 2.5D, ideal 
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Figure 15. The change in total energy of the system is plotted 
against time for hall-med-hr (thick line, upper panel) and hall- 
high-hr cases (thick line, lower panel). In each case the change in 
total energy is plotted for fuU-low-hr (thin line) for comparison 
and also for hd-zoro-hr (dashed line). 
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Figure 16. The total perturbed magnetic energy in the system 
is plotted against time for fuU-low-hr (thin line) and hall-high-hr 
(thick line) cases. The magnetic energy is seen to experience very 
different growth in hall-high-hr. 



MHD system is made possible b y the Hall ter m introducing 
a handedness into the flow (e.g. IWardlelll999l ). 

If we follow the dynamics further we find that the KH 
instability, in the presence of high Hall resistivity does not 
saturate to a quasi-steady state as it does in, for example, 
the fuU-low-hr case. As the z components of the current and 
magnetic field continue to grow, the Hall effect now acts 
on the non-parallel currents and magnetic fields that have 
arisen between the z-directions and the sy-plane. This has 
the result of re-orienting some of the magnetic field, and 
electron fiuid flow, back onto the xy-p\a,ne. During this pro- 
cess the electron fluid obtains a velocity away from the KH 
vortex, which results in a broader volume of plasma being 
disturbed. This feeds the continuous growth of the magnetic 
energy in the xy-plane, and thus causes continuous growth 
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Figure 17. The perturbed magnetic energy in the system is plot- 
tod against time for hall-high- hr (solid line). Also plotted for com- 
parison are the growth of the magnetic energies in the xy-plane 
(dashed line) and in the z-direction (dot-dashed line). It can be 
seen that there is significant growth in the z-direction. 
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Figure 18. Plot of the magnitude and vector field of the magnetic 
field in the xy-plane in hall-high- hr at time 8 4s. It can be seen 
that the magnetic field does not undergo as much wind-up as is 
seen in the bulk velocity field (upper panel of figure [TOt . 
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of the electron transverse kinetic energy, as can be seen in 
figure [20l In this particular simulation, the size of the simu- 
lation grid is set so that it provides a sufficiently large region 
of ordered flow that can be transformed into disordered flow 
during the time of the simulation. The instability will nat- 
urally saturate when it has exhausted the area of ordered 
flow available to it. 

To summarise, through their strong decoupling from the 
magnetic field, the dynamics of the bulk fluid and ion fluid 
demonstrate behaviour very similar to that of hd-zero-hr in 
which the KH vortex remains intact. The high Hall case 
can in fact be thought of as two separate systems occur- 
ring simultaneously; the bulk fluid demonstrating hydrody- 
namic behaviour and the continuously widening volume of 
perturbed fluid perpetually feeding the growth of the mag- 
netic field through the Hall effect. These two systems are in- 
trinsically entwined through the requirement of charge neu- 
trality, by which the electron fluid causes a widened area 
of perturbed ion fluid, and thus the bulk fluid. Both sys- 
tems are relatively energy efficient and, following the initial 
growth of the KH instability, the overall system experiences 
little further loss of energy. The supply of energy to feed 
the magnetic field is limited only to the physical size of the 
computational domain over which the simulation is run. 

This dynamo action occurs only under certain condi- 
tions. As the initial Hall resistivity is increased from moder- 



Figure 19. Plot of the magnitude and vector field of the bulk 
velocity (upper panel), ion velocity (middle panel) and electron 
field (lower panel) for hall-high-hr. The dynamics of the electron 
fluid can be seen to be very different to the other two fluids, and 
more similar to the magnetic field (figure (TSj. 




Figure 20. The transverse kinetic energy of the electron fluid 
is plotted against time for fuU-low-hr (thin line) and hall-high-hr 
(thick line). Also plotted in the electron transverse kinetic energy 
for hd-zero-hr (dashed line). 
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ate to high, the electron fluid, and thus the magnetic fleld, 
becomes increasingly decoupled from the neutrals. As a re- 
sult, the neutral fluid tends toward hydrodynamic behavior. 
This dynamo action is observed only when the Hall resis- 
tivity is increased to a sufficiently high value that the KH 
vortex formed by the neutral fluid is no longer constrained 
by the magnetic field. Unlike the pure hydrodynamic case, 
the presence of charged fluids undergoing different dynamics 
leads to the dynamo action observed. If the Hall resistivity is 
not sufficiently large, the magnetic field eventually leads to 
the destruction of the KH vortex in the neutral fluid through 
reconnection as in the ideal MHD case. 

Even a moderate amount of Hall resistivity results in 
a wider volume of fluid being disturbed by the instability. 
This agrees with pr evious studi es of the Hall effect on the 
KH instability (e.g. lHuballl994 'l. The re-orientation of the 
magnetic fleld lines within the KH instability has al so been 
obser v ed in s tudie s of MRI in accretion disks (e.g. IWardld 
I1999I ) . iKunj l|2008l ) investigated a simple model of accretion 
disks without rotation and demonstrated that the combined 
actions of the shear instability and the Hall effect leads to 
increased stretching of the magnetic field lines. This was 
shown to result in continued growth of the instability, which 
corresponds well to the dynamo action observed in our sim- 
ulations with high Hall resi stivity . It is imp ortant to note , 
though, that studies by both lKunj l|2008h and lWardld l|l999l l 
are linear studies and, as such, do not extend to the nonlin- 
ear r«girnfi_of_theinstability. A more recent numerical study 
bv lNvkvri fc Ottd (|200J) demonstrates the twisting of mag- 
netic field lines, but due to the inclusion of magnetic re- 
connection, doesn't produce the magnetic dynamo observed 
here. 



5 CONCLUSIONS 

We have presented the results of a suite of fully multifluid 
MHD simulations of the KH instability in weakly ionised fiu- 
ids such as, for example, molecular cloud material. Through 
varying the collision coefficients between the various charged 
species and the neutrals we were able to investigate systems 
in which ambipolar diffusion dominates the multifluid effects 
and ones in which the Hall effect dominates. We validated 
our KH simulations through comparison of an ideal MHD 
simulation with previously published results and performed 
resolution studies for each of these cases to ensure that our 
conclusions are not unduly effected by our numerical reso- 
lution. 

Our findings can be summarised as follows: 

• The multifiuid effects do not significantly influence the 
linear growth rate of the instability. 

• Ambipolar diffusion dramatically reduces the energy as- 
sociated with the perturbed magnetic fleld. This happens 
through diffusion for moderate ambipolar resistivity, but 
through both diffusion and decoupling of the magnetic fleld 
from the bulk flow for higher resistivity. 

• The Hall effect, as expected, rotates the magnetic fleld 
out of the initial xy-p\a.ne. 

• In contrast to both ambipolar dominated and ideal 
MHD flows, the Hall effect causes the system to fail to settle 
to a quasi-steady state after saturation of the instability. 



• For moderate Hall resistivity the perturbed magnetic 
fleld contains higher energy than in the ideal MHD case, 
presumably due to a fleld topology which impedes (numeri- 
cal) reconnection. 

• For high Hall resistivity strong dynamo action is seen 
as energy associated with the magnetic field grows without 
any apparent signs of saturation. 
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